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1 Abstract 

We study the first steps of ice formation in fresh water turbulent convection (frazil ice regime). 
We explore the sensitivity of the ice formation process to the set of non-dimensional parameters 
governing the system. We model the mixture of ice crystals and water as a two-phase medium 
composed of water and ice particles of fixed diameter. We use the Boussinesq approximation and 
we integrate numerically the set of equations making use of a numerical code based on second 
order finite difference. A dynamic LES model for the subgrid scales is used. 

We show that the ice particle rise velocity and the ice concentration source term coefficient 
significantly affect the frazil ice dynamics. The maximum of ice production is obtained in those 
situations where the rise velocity is of the same order of magnitude of a characteristic velocity 
of the thermal downwelling plumes. We develop a simple model which captures the trend of the 
growth rate as a function of the relevant parameters. Finally we explore the nonlinear regime 
and we show that the parameter which plays a key role in fixing the concentration value at the 
statistically steady state is the heat exchange source term coefficient. 

2 Introduction 

The first stage of ice formation at the supercooled free surface of oceans, rivers and lakes is 
called frazil ice, a suspension of individual, small randomly-oriented crystals typically measuring 
approximatively 1-^4 mm in diameter and 1 100 /iin in thickness |10 | .|13 | . 

Frazil ice evolves rapidly into a thin slurry of ice platelets giving the sea a greasy or milky 
appearance, known as grease ice. As cooling continues and the density of the frazil suspension 
increases, the frazil crystals freeze together into small cakes. Collisions between neighbouring cakes 
forces frazil ice over the cake edges where rims form to give the characteristic form of pancake, 
hence the appellation pancake ice. Initially formed pancakes may be only a few centimeters in 
diameter, but they grow in diameter and thickness from surrounding frazil matrix and may reach 

3 -7- 5 m in diameter and 50 -f- 70 cm in thickness. The last stage of sea ice evolution is the 
coalescence of pancake ice floes, which aggregate into a continuous ice sheet. In figure QJ a 
picture of a frazil-pancake ice field in Polar ocean is reported. 

In rivers, frazil ice production is a transient phenomena which occurs in turbulent supercooled 
water. The first stage is the water temperature lowering to about 0.01 0.1 °C below the freezing 
point because of the wind action or the radiative cooling; then the frazil discs start to form and 
the water temperature following the mass of crystals increases in the order of minutes toward the 
freezing point. The presence of frazil ice can cause serious problems to hydroelectric facilities such 
as the blocking of turbine intakes, the blockage of hydroelectric reservoirs and the freezing open 
of gates. 

In polar ocean, frazil ice can form under ice covers and in open region within polar pack ice, 
more precisely in 100 meters scale leads and in kilometer scale polynyas. Frazil ice is typically 
present in the Marginal Ice Zone (MIZ), which is the transition region between the open polar 
ocean and the continuous ice that covers the central basin. Due to the combined action of wind 



and surface waves, the turbulent waters of MIZ prevent the formation of large ice sheets. The 
presence of frazil ice fields is associated with important geophysical and biological processes of 
the polar oceans so that the study of frazil ice dynamics is an active research field. As frazil ice 
production is accompanied by salt rejection, it is believed that it could play an important role in 
stimulating the convection process of ocean waters. Moreover, frazil ice plays an important role 
as scouring and sediment transfer agent and also it can transport biological species from the sea 
bottom to the surface pack ice. 

This paper is organized as follows. Section is devoted to a brief description of the frazil ice 
formation process. In section 0] the derivation of the frazil ice transport model is reported, while 
section [5] is dedicated to the Boussinesq-like approximation. In section we present numerical 
results of the simulations and finally section [7| is devoted to the conclusions. 




Figure 1: Frazil-pancake ice field in Polar ocean (www.vims.edu/bio/microbial/NBPice.html). 



3 The formation of frazil ice 

The time it takes for ice formation in oceans, rivers and lakes can be quite different from one year 
to the other, due to meteorological variations. Estimating this time is a very complicated task, 
which involves many parameter such as wind, water temperature and water depth. Immediate 
applications of this problem can be found in shipping and ice breaking service. To this end, 
many theoretical and experimental studies were performed during the last decades, giving useful 
information about time histories of temperature and ice formation. 

Designers of hydraulics structures such as canals, water transport facilities and hydroelectric 
power structures are particularly interested in three ice-flow regimes, which are common in rivers 
during freeze-up period: continuous ice sheet, floating frazil ice and well mixed flow (with frazil 
ice crystals vertically mixed in the water flow). In the last case, if the velocity of the river is 
sufficiently large, open water conditions may persist throughout the winter. In consequence of 
this, a very important task is to predict, starting from a given set of flow conditions, which kind 
of regime could be expected, and the transition from an ice cover to open water. Obviously, when 
frazil ice crystals remain well mixed with the fluid, they cannot accumulate beneath the surface 
and the ice cover can not form, so the conditions under which the transition from well mixed 
flow to floating ice layer occurs, could be employed as a rough criteria for the formation of the 
continuous ice cover. In the rough criteria cited above was compared, with good agreement, 
with experimental and field data. In particular, it has been remarked that the presence or absence 
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of the well mixed regime could be forecasted comparing the buoyancy time scale Tg with the 
vertical turbulent diffusive mixing scale Tjj. When Tg < To a layered flow will develop, because 
the frazil ice crystals rise to the water surface faster than they are removed by turbulent mixing. 
Under the same argument, the opposite case Tg > Tr> corresponds to a well mixed flow. More 
sophisticated criteria take into account other hydrological aspects, such as sinuosity, channel and 
slope variability, the transition from frazil ice layer to ice cover throughout pancake ice formation 
and the mechanics of interaction between the different kinds if ice [2] ■ 

In |15j a numerical simulation of ice formation and transport in turbulent ocean water was 
performed. The upper ocean was modelled as an Ekman layer, using a k — e turbulence model 
for the Reynolds-averaged Navier-Stokes (RANS) equations. This model was improved in [T§| 
considering the interaction between frazil ice particles of different size. In [§] a model of frazil ice 
transport under an ice shelf was proposed and a Boussinesq-like approximation was adopted in the 
limiting case of very small ice concentration. In this model the vertical structure was simplified 
by assuming all properties to be constant over the thickness of the plume and discontinuous at its 
boundaries; the problem was reduced to one dimension. In |18| the model was extended to the 
case of ice particles of different size. More recently |TJ|, a Large Eddy Simulation (LES) with a 
Smagorinsky turbulence model was carried out in order to study frazil ice dynamics and the ice 
effects on turbulent structures under polar ice covers and in leads. This work took into account 
the effects of salinity. Finally |Hj presented results of LES simulations for a turbulent flow with 
frazil ice crystals of different size. 

In our work we propose a parametric study of the problem in the framework of Boussinesq 
approximation. In the case of large ice concentrations, the particles interaction at the ice crystal 
microscale should be considered. In particular, it was observed that this phenomenon becomes 
relevant for ice volume concentrations above .003 |15j . Moreover, the presence of frazil ice crystals 
produces an increment in the ice- water mixture effective viscosity. In [3] it was shown that, for 
small ice concentrations, this parameter differs from that of pure fluid by a value of the same order 
of the concentration itself. 

In this paper we have used a convection model for frazil ice transport in turbulent water by 
considering the mixture of ice crystals and water as a two-phase medium. The balance equations 
for mass, momentum and energy of the ice-water mixture are supplemented by an equation for the 
temporal evolution of the ice mass concentration. This equation is strongly coupled to the others. 
In particular, the ice particles melting/freezing process and the temperature fluctuations modify 
the density of the mixture and thus the buoyancy force, which, in turn, affects the hydrodynamical 
behavior. 

The simulation of these equations is subject to numerical instability due to large round-off 
errors. In order to avoid this behavior, related to the large magnitude of the coefficient of the 
buoyancy term, we have applied a Boussinesq-like approximation to the equations of motion. 

A complete model of frazil ice dynamics should be able to represent the complex processes 
involved in the ice formation, like growth of ice crystals, flocculation (the clustering of two, or 
more, frazil particles), break-up (the cracking of an ice particle) and secondary nucleation (the 
production of small ice crystals by the ice fragments scavenging from crystal surface). As a 
consequence, a large range of particles size distribution is expected with the crystal diameters 
varying from micrometer to millimeter. In this paper, however, we have considered only one 
particle size. This is not a crude simplification, because the small ice crystals rapidly evolve into 
larger ones and finally only ice particles within a narrow size range actively affect the system 
dynamics HJ. 

For the sake of simplicity we limit our analysis to ice formation in fresh water and we postpone 
studying the effect of salinity to a future work. 

We have integrated numerically the evolution equations in the Boussinesq limit using a second 
order finite difference method. The smallest scales have been modelled within the Large Eddy 
Simulation approach: the equations for the resolved scales are obtained introducing a spatial filter 
that removes the unresolved small scales of turbulence; the filtered nonlinear terms are modelled 
using the approach of |20|. Numerical simulations have been carried on to determine the sensitivity 
of ice growth rate and saturation process to different sets of non-dimensional parameters. 
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4 Model equations for frazil ice transport in turbulent wa- 
ter 



We model the convective motion in the upper ocean mixed-layer as convection between two parallel 
and infinitely long boundaries (albeit periodic), with the upper boundary corresponding to the 
air-ocean interface and the lower one corresponding to the mixed layer bottom. The distance 
L between the boundaries and the temperature difference AT maintained between them are set 
equal, respectively, to a typical depth of the mixed-layer (30-^100 m.) and to a typical temperature 
gap between the air-ocean interface and the mixed-layer bottom (0.1 -j- 0.2 °C). The temperature 
at the air-ocean interface is set to a value slightly less than the salt water freezing point — 2.0 °C 
corresponding to a typical salinity value So = 34.5 psu. Salinity variations play only a small role in 
the thermodynamics of the system and we shall hold it constant [H] . Ice concentration is assumed 
to be large enough to be represented as a continuum but small enough to neglect the particles 
interactions. As a consequence, the system of water and frazil ice particles can be modelled as a 
two-component mixture, which is treated as a locally homogeneous fluid with averaged properties. 
In the dilute limit that we are considering, the effective viscosity of the ice-water mixture differs 
from that of pure fluid by a value of the same order of the concentration itself and this viscosity 
increment can be neglected. 

We consider a frame of reference such that the z axis is vertically upward. The ice-water 
mixture density p is given by: 



P = Pi<fii + Pw^w, (1) 

where pi and <pi are, respectively, the ice density and volume fraction; similarly p w and <f> w = 1 — 0, 
are, respectively, the water density and volume fraction; p w is related to the temperature by the 
following equation of state: 



Pw = Pwo[l - a{T - T )], (2) 

where p w o is the density corresponding to the lower boundary temperature To and a is the thermal 
expansion coefficient. 

The total density p can be written as: 



V pi p,„ I 



in terms of Ci, that is the ice mass concentration, obeying Cip = pi4>i (similarly C w = 1 — Cj will 
denote the water mass concentration). 

The equation of continuity for the total density p is given by: 

^ + V-(pu) = 0, (4) 
while the equation for the mixture momentum balance is formulated as follows: 

p-^+p(u- V)u= -Vp + pgk + vpV 2 u, (5) 

where u and p are respectively the velocity and the pressure fields, g is the gravity acceleration 
and v is the kinematic viscosity. 

The equation for the temperature T can be written as follows: 

dT 

— + u-VT = fc T V 2 T + G T , (6) 

where kx is the water thermal diffusivity and Gt is the source term due to melting and freezing 
processes, considering the ice particle as a sphere: 
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AT 

G T = — = X T pC i (T i -T), (7) 



where 7j : 

3h 



At = 



RsCpPwPi 

and h = is the heat transfer coefficient (Nu is the Nusselt number and k w is the sea water 

thermal conductivity, R$ is the ice particle radius), Tj is the ice particle temperature, taken equal 
to the freezing point, and c p is the water specific heat. 

Finally, the equation for the ice component can be written as ([5],E2> E3) : 

+ V • (upG) = w r ^^- + fc 4 V • (pVC 4 ) + d, (8) 

where ki is the ice molecular diffusivity, w r = (0, 0, w r ) is the rise velocity that derive from the 
balance between the upward net buoyancy force and the drag force. It is still an open question how 
to model this velocity; we have left it as a free parameter for sensitivity analysis f[T].|5].[T5 ] .|19 j ). 
Here G,; is the source term due to melting and freezing processes and it is given by: 

G i = ^p-=\ iP C i (T i -T) (9) 



where [7]: 



A; 



3ft 



RsLiPi 

and Li is the ice latent heat of fusion. 

It is important to remark that we have used spherical particles instead of disk-like ones which 
better represent the frazil ice platelets. The choice is motivated by the great simplifications 
introduced in the model. Moreover the ice particle shape affects the source terms only in their 
magnitude; the latter parameters will be used in the foregoing for the sensitivity analysis and they 
will vary in a wide range; so we can conclude that the effect of the particle shape on the source 
terms is not significant for the analysis we are going to carry on. 

In order to rewrite eqs. (1211 )[l in non-dimensional form, we consider the depth of the upper 
ocean mixed-layer L = L as length scale, so the depth of the layer is 1, the viscous scale velocity 
U = vjL as velocity scale, the corresponding time t — L/U — L 2 /v as time scale, the water 
density p — p w o as density scale and the temperature gap AT between upper and lower boundary 
as temperature scale. We indicate by z the dimcnsionless vertical coordinate with zq = —1/2 the 
lower boundary. We obtain: 

p w = l-a(T- T ), (10) 



P : 



V Pi Pw ' 



^ + V-(pu) = 0, (12) 

9u , _ Ra -> ~ . . 

p— + p(u ■ V)u = -Vp - — pk + P V 2 u, (13) 

dT 1 ~ 

— +u- VT= — V 2 T + /3 P a(T-T), (14) 

^1 + V • {upd) = + J^V • (pVQ) + 1P C l {% - T), (15) 
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where Ra = gaL 3 /(kv) is the Rayleigh number, Pr = u/k is the Prandtl number, Sc — vjki 
is the ice Schmidt number, w r = w r /U, a = a AT, = ^f 1 , 7 = A '^ r ' , T = T/AT and 
fb = Tq/AT. 

At the lateral boundaries, we prescribe periodic conditions, while at the horizontal boundaries 
we impose a fixed temperature and free-slip condition for velocity. For frazil ice concentration we 
impose Dirichlet conditions Cj = at the bottom boundary, while the conditions to impose at the 
top is still an open question; we postpone the discussion of such boundary conditions at section 
IO and IOI 



5 The Boussinesq approximation for frazil ice transport 
model 

The buoyancy term in the mean momentum balance equation (|13J) is multiplied by the factor: 

Ra 



aPr 



In the cases of physical interest, the Rayleigh number is in excess of 10 8 -7- 10 9 and a ~ 10 -6 ; 
hence this factor can be estimated to have a magnitude of at least 10 14 -j- 10 15 . As a consequence, 
the numerical evaluation of this term in a double precision code is subject to significant round-off 
errors: this produces severe numerical instability. In order to avoid this problem we introduce a 
Boussinesq-like approximation. We introduce a hydrostatic equilibrium without ice and then we 
write down the equations for the small perturbations to this state. The fluctuations are written as a 
power series of the small parameter a and the expansion is truncated to lowest order. Temperature 
fluctuations, pressure and velocity field are of order 0(1), whereas density and ice concentration 
fluctuations are of order 0(a). We thus have: 



+ 0(a), Vp 



Ra 
aPr 



[I + a(z - z )]k + Vp' + 0(a), 



T = T o -(z-z o )+T' + 0(a), Q = aC[ + 0(a 2 ), p = 1 + a(z - z ) + ap' + 0(a 2 ). 



Substituting into eq. I|ll|l to leading order in a we obtain: 

l -^ic' i + T' 
\ Pi 

and the equations (|12J) and Ijl3(l become: 



p' = 



V • u' = 0, 



(16) 



(17) 



^ + (u'-V)u' = -Vp'+- 



1 - P i 



C' + T' 



VV. 



(18) 



Note that the factor -M^ in the buoyancy term of eq. I|13fl has been transformed to , which is 
not subject anymore to numerical instability due to round-off errors. 
Temperature equation (jl4[) becomes: 



BT' 1 

^+u'-VT' = -w 1 + — V 2 T' + fiaCjiz, T% 



(19) 



where: 



f(z,T')=T i -T -(z-z )-T' 
and finally the ice transport equation l|15|) reduces to: 



G 



^ + u' • VCI = jr^C* + w r ^ + 7 C'f(z, T'). (20) 

Equations (|17I20[) are a complete set of equations for the small perturbations to the hydrostatic 
equilibrium. For sake of simpler notations, we shall drop from now on the primes from the variable 
names. 



6 Numerical results 

Equations I|17I20[) were solved numerically using a fractional step approach |12j with a third 
order, three steps, explicit Runge-Kutta scheme. For the spatial discretization, we have employed 
a second order accurate finite difference scheme on a staggered grid. For the momentum equation, 
each Runge-Kutta time step has been splitted into two steps; the first one consists in the calculation 
of a temporary velocity field without making use of the pressure terms; the second step is a 
projection of the solution onto a divergence-free space, which results in solving a Poisson equation 
for pressure. A fast Poisson solver using Fourier and cycling reduction method was adopted |16| . 

The effects of the small unresolved scales on the resolved ones are modelled by using the 
Subgrid Scale Model proposed by [20! ■ We utilize the same subgrid scale model for the flux of ice 
concentration and for the temperature: the turbulent Schmidt number is computed by using the 
dynamic procedure and varies in space and time. 

The domain has an extension of 27r in the horizontal directions and height equal to one. The 
resolution varies with Rayleigh number: at Ra = 10 4 , 10 6 and 10 8 we used respectively 31 x 31 x 30, 
51 x 51 x 36 and 81 x 81 x 60 collocation points. 

6.1 Ice production in deep water 

A significant amount of oceanic frazil ice might be produced in the deep layers of water. In |13| . 
for example, it is stated that a large amount of frazil ice in seawater adjacent to shelves could be 
produced by freezing of parcels of water rising towards supercooled upper layers. In this series 
of simulations, we want to focus on the ice produced in the regions far from the upper layer. To 
this end we use the boundary conditions Ci = at z — ±1/2. These conditions amount to cancel 
the convective flux of ice through the horizontal boundaries. They produce a very thin boundary 
layer immediately beneath the free surface where a non-physical ice concentration drop is present. 

In all the subsequent simulations we assume that the freezing temperature is T; t = T + z Ql that 
is the upper half of the hydrostatic layer is below freezing point. 

Typical curves of the total ice concentration versus time are shown in Fig. |3 In the left figure, 
we clearly distinguish a short transient followed by an exponential growth, ended by the nonlinear 
saturation. Evolution of kinetic energy is shown in Fig. The ice seeds are introduced at time 
t = 0. We observe that the production of ice is fast compared to the large eddies turnover time. 
As a consequence, during the linear phase, the velocity and temperature fields are almost constant 
and we get an almost pure exponential growth of ice concentration. In the right figure, instead, 
the exponential growth phase does not appear very clearly. This is because the time scale for the 
ice formation is large compared to the convection turnover time: the flow changes significantly 
during ice formation. The fast and slow regimes appear very clearly in Fig. 

In the fast growth regime, we observe that most of the ice is produced inside the strongest 
downwelling plumes. A typical situation is shown in Fig. it is a snapshot of the flow taken during 
the linear phase of ice growth. The isovalue for temperature fluctuation has been chosen so that 
only the most intense downwelling plume appears (since the flow is periodic in horizontal directions 
the four large green structures are actually contiguous). We see that ice (white structures) is 
present mostly at the top of that plume. This is natural, since it is also the coldest place inside 
the flow: the largest part of ice is produced there, and, since the time scale for the linear phase 
is small, the shape of the plume almost does not change and the ice remains confined inside it. 
By contrast, in the slow regime, the spots of large ice concentration do not necessarily reside in 
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Figure 3: Evolution of kinetic energy for the two simulations of Fig. Ice seeds are introduced 
at time t = 0. Full line: 7 = 10 4 , w r — 150 (fast growth). Broken line: 7 = 10 3 , w r = 10 (slow 
growth). 
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the large downwelling plumes. Fig. [S] shows a typical situation where we see small spots of ice 
next to the strongest downwelling plumes, and a large spot in a relatively quiescent region in the 
flow. The ice still forms mostly inside the coldest downward plumes, but now the process of ice 
formation is slow in comparison to the convection time scales and the ice forming has time to 
migrate towards quiescent regions. Looking again at Fig. |3| we remark that the total kinetic 
energy decreases as soon as a significant amount of ice is produced; this happens because most 
of the kinetic energy is concentrated in the most intense plumes: the ice that is produced in the 
strong downwelling plumes modifies the buoyancy force and tends to slow down the downward 
motion; as a consequence, the kinetic energy of those plumes is reduced. 



Figure 4: Isosurface of negative temperature fluctuations (-0.2, green) and isosurfaces of large 
ice concentration (.lC max , white). Snapshot taken at a temporal station during the exponential 
growth regime for a fast ice formation regime (Ra ~ 10 6 , 7 = 10 4 , w r — 150). 




In Figs. HO and we plot the growth rate for the linear phase as a function respectively of the 
rise velocity w r and of the thermal exchange coefficient 7. Each point in the figures is obtained 
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by carrying on a full simulation until nonlinear saturation of ice concentration occurs. In the fast 
regime, the flow is almost frozen during the ice production and the growth is almost exponential; 
the growth rate is computed as the slope of ice concentration versus time in log scale (as in Fig. O 
(left)), without considering the initial transient and the final nonlinear regime. In the slow regime 
the value of the growth rate is determined as the slope of the line which fits at best the curve in 
the linear phase, omitting the initial transient. 

We notice in Fig. H3 that the growth rate has a maximum at a particular value of the rise 
velocity which we denote W. For Ra = 10 6 it is about w ~ 75. 

Looking at the concentration profiles averaged in the horizontal directions (Fig. [SJ we can 
distinguish three situations according to the value of w r . 

In the large rise velocity regime (w r > w) , the ice crystals are rapidly pushed upward and they 
cannot act as nuclei for ice formation; therefore the ice rate of growth decreases when rise velocity 
increases and most of the ice concentration is found near the top. If w r ~ w, the rise velocity is 
balanced by the downwelling velocity of the plumes in the upper layer of the domain, where the 
temperature is the lowest. Crystals remain confined in this region and their concentration grows 
rapidly. Finally if w r < w the ice is transported down and exits the plumes where it was created. 
The ice spreads throughout the domain and the average concentration profile is almost flat (Fig. 

03). 

Note that the profiles observed by ([21]) beneath lead in the Beaufort Sea (north of Alaska) and 
those obtained by f |17p from the numerical simulations of lead dynamics show a large presence 
of ice near the top of the layer; this is in agreement with the situation where w r > w. 




rise velocity (w r ) coefficient y 



Figure 6: Growth rate of ice concentration Figure 7: Growth rate of ice concentration 

versus rise velocity for Ra = 10 6 and 7 = versus 7 for Ra = 10 6 and w r = 150. 

10 4 . 



6.2 A ID model for fast ice formation 

We propose in this subsection a simple model that explains the behavior of the curves shown in 
Figs. El and for the fast regime. We start from {2D)) and assume that inside the intense vertical 
plumes we can neglect horizontal velocities and horizontal derivatives. Therefore Cj is a function 
of z and t and we get: 

W = w + w r , r=7(2i-T). (22) 



where: 

1 1 1 



S Sc Sct ur b 

We perform an analysis into normal modes: 



Ci(t,z)=C(z)e u 
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Figure 8: Average ice concentration for different values of w r (Ra — 10 ,7 = 10 ) 



where u) € C is the eigenvalue and C(z) the eigenvector: 

C"{z) - WSd'(z) + (r - uj)SC{z) = 



(23) 



(C and C" denote derivatives with respect to z). Let us further assume that W, T and S 
are constant. We look for solutions with exponential dependence on z. Imposing the boundary 
conditions: 



C(-l/2) = C(l/2) = 0, 



we get: 



uj = r- 



W 2 S 



S 



|n| ^ 



Ci(z,i) = ae^e^ ([1 - (-1)"] cosrmz + i[l + (-1)"] sinmrz) , 

where n is a non vanishing integer and a is an arbitrary constant. The case n = is degenerate; 
the growth rate is: 

W 2 S 

u=T-^f-. (24) 
We can now impose only one boundary condition. If W > the eigenvector is: 

" >s WSz 

C{z) =ae—(l-2z) 

and the boundary condition is exact in z = 1/2 and only approximately zero in z = — 1/2. If 
W < 0: 

C(z) = ae^(l + 2z) 

and the boundary condition is exact in z — —1/2 and only approximately zero in z = 1/2. 

The degenerate case is the one with the largest growth rate and it is therefore the most 
interesting. In Figs. HOI and lTTI we plot the growth rate to versus respectively w r and Y together 
with the actual values obtained by the LES simulations. The values of W, T and S are expressed 
by (1221) and depend on typical values of vertical velocity w, temperature T and turbulent Schmidt 
number Scturb- Those values could be predicted by turbulent convection theories but we decide 
here to choose the values which fit at best the LES curve. As we see in Fig. ^| the fit is very 
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good for rise velocities greater or equal to w. The bad behaviour at smaller values is due to the 
fact that thei ice does not stay confined inside downwelling plumes and therefore the hypotheis of 
our model is not verified. In order to get a more quantitative agreement, we refine the model by 
dropping the assumption that W and T are constant. We compute the actual vertical profile of 
these functions by averaging the values obtained by the LES simulations inside the most intense 
downwelling plume (see Figs. EJ. 




(a) (b) 

Figure 9: Vertical profile of temperature fluctuations (a) and vertical velocity (b) of the two most 
most intense downwelling plume. Data are averaged over the horizontal dimension of the plume. 
Ra = 10 6 , 7 = 10 4 , w r = 150. 

The differential eigenvalue problem becomes: 

-=C" - W(z)C' + T{z)C = cud, 

(25) 

C(-l/2) = (7(1/2) = 0, 

for prescribed functions W{z) and T(z). The solution is no more analytical and we compute 
it numerically. We adopt a spectral technique: we expand the solution in series of truncated 
Chebyshev polynomials and we discretize the problem over the Gauss-Lobatto grid points: 

1 7TJ 

Zj = - cos — j = 0, . . . , N. 
1 2 N 

We obtain a linear algebraic generalized eigenvalue problem: 

AC = ujBd, 

where A and B are (N+l) x (N + l) matrices. A method to solve such problem is to use 
the QZ algorithm. Such algorithm gives the full spectrum of eigenvalues/eigenvectors. Another 
method is to use iterative techniques (Krylov based methods), such as the Arnoldi method or the 
unsymmetric Lanczos method 0- This method gives few eigenvalues corresponding to the least 
stable part of the spectrum. For our problem we have implemented both the QZ and the Arnoldi 
method and we have checked that the results produced by both techniques are always similar. 

The growth rates obtained by solving H25|) with the vertical profiles of the two most intense 
downwelling plumes and those given by the constant coefficient model l|24|) are shown in Figs. 
ITU1 and ^2 From Fig. ^] we observe that when the rise velocity is far from the peak value W, 
the growth rate predicted by our model underestimates the actual value (this occurs for about 
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w r > 300 and w r < 75). The reason is that when the growth rate is small the flow cannot be 
considered steady during the exponential ice growth, but rather it fluctuates. 

Figures IT21 and IT31 are similar to Figs. 1101 and lTTl except for the Rayleigh number that is 10 8 
instead of 10 6 . We see that the results obtained at Ra = 10 6 are still valid, thus they are robust 
and we can hope to extrapolate them to even larger Rayleigh numbers that are of interest for 
geophysical applications. 

Some analytical progress far from the fast regime is possible provided w r is not too small (so 
to be able to balance the downwelling in the plumes). In this way, the permanence time t p of an 
ice crystal in a plume can be considered long and the growth process will be dominated by the 
plumes with the largest growth rate: (exp(rr p )) ~ exp(r MAX r p ). Taking at the same time the limit 
Tt <C 1 with t the typical plume circulation time, allows to treat the velocity and temperature 
fluctuations as a noise and to calculate averages in closed form. 

Let us look at a single Lagrangian parcel, moving with velocity w(z(t), t) + w r (plus the effect 
of molecular diffusion) and endowed with an ice concentration c(t) obeying: 

^- =1 [ Ti -T(z(t),tMt). 

Identify with tilde and overbar, the fluctuating and mean field component of w and T at the parcel 
position z: w{z, t) = w(z, t) + w(t); T(z, t) = T(z, t) + T(t) In the limit r — > 0, the equations for 
z(t) and c(t) can be written in the form: 

dz = [w r + w(z,t)]dt + £7,jT 1/2 dBi + (2/Sc) 1 / 2 dB 2 , 
dc = 7[Tj - T(z,t)]cdt + -yafT^cd^Ba, 

where the dBi are Brownian increments, that, for simplicity, we take as independent: {dBi) = 0, 
(dBidBj) = Sijdt. The superscript in d (S) i33 indicates Stratonovich prescription 0], and arises 
from the short correlation time limit f* +T c(t')f(t')dt' ~ icr f .T 1 / 2 [c(t)+c(t + r)][B 3 (t)-B 3 (t + r)]. 
This leads to a correction to the mean growth rate, and the second of eq. (|6.2|) will take the form, 
back to the standard Ito prescription: 

dc = 7 [Ti - f(z, t) + \pa\r\cdt + ja t T^ 2 cdB 3 . 

The PDF p(z, c; t) will then obey the Fokker-Planck equation: 

r)n r) 8 1 fl 2 1 Ft 2 

P + J Wr + u)p} + 7 -[(T, -T+ -^ 2 f r)cp] = -(p/S) + --( 7 4r C 2 p), (26) 



where 5 _1 = Scq 1 + o\t jl may include non-subgrid contributions to the fluctuations. The 
concentration field resulting from this averaging procedure will be: Ci(z, i) — J p[c, z; t) ede, and, 
substituting into cq. 12(j|) . we obtain: 

FlC rl f) 2 - 1 

^ + £[K + w)Ci] = ^Ci/S + 7 (T ( —f+ -l4r)a. (27) 
Thus, the temperature fluctuations lead to a net increase in the ice concentration growth rate: 



r-r+- 7 4r. (28) 



Neglecting this effect, as shown in Fig. leads to underestimating the actual growth rate. 
Similarly, in Fig. II II the slope of theoretical curves should be increased, which is precisely what 
eq. (J2E1) does. 

Notice that Eq. (|21l) would be recovered from 1271) (as it should) if no fluctuating component 
were singled out of the w and T fields. 
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rise velocity (w ) 



Figure 10: Growth rate of ice concentration versus rise velocity for 7 = 10 4 and Ra = 10 6 . 
Full line: Growth rate from LES simulations. 
Broken line: eigenvalue given by l|24|) . 

Line with circles: eigenvalue of problem l|25(l for the most intense downwelling plume. 

Line with triangles: eigenvalue of problem i|25|) for the second most intense downwelling plume 
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Figure 12: Growth rate of ice concentration versus rise velocity for 7 = 10 5 and Ra = 10 : 



~ le+05 
5 



50000 









1 1 


1 ' 




LES 


y 




eigenvalue 






— plume (1) 


s . 
s 




a — a plume (2) 


X _/ y 
* y$ 
✓ y^^.y 






S //*' 
^ S^oS y 






* /V,-' -fir 






S />- 
✓ >yyS -fir 






* y£r^ 
y sgf ^ 

y 

✓ s^*y^ 
* yr -y 
y r> ^^r 

y -^r^ 






* >y^y^ 
' y^ 

y — 

l.l.l, 



le+05 2e+05 3e+05 4e+05 5e+05 

coefficient y 



Figure 13: Growth rate of ice concentration versus 7 for w r = 1000 and Ra = 10 8 
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6.3 Non linear phase 



The nonlinear saturation occurs when the ice concentration (rescaled by the factor a, see ((201) 
reaches values of the order of unity: the density of the mixture varies by amounts of order a and 
the buoyancy force in 1|18() is significantly affected. One must bear in mind that the Boussinesq 
approximation breaks down when ice concentration exceeds values of order unity; therefore results 
of the nonlinear stage must be taken with caution: they give qualitative indications but they 
cannot be used to assert quantitative conclusions. 

In Fig. 1141 we plot the temporal evolution of total ice concentration for three simulations with 
different values of (3. The linear stage is insensitive to the value of (3 because the ice concentration 
is too low to produce any change in temperature fluctuations. However we see that the nonlinear 
saturation is affected by the value of (3 and the general trend is that for larger values of j3 the ice 
saturation levels tend to decrease. 

In Fig. E| we see that the presence of ice inhibits convection. 

In Fig. 1161 we plot the average vertical temperature profiles during the nonlinear stage. We 
see that temperature tends to increase when (3 decreases. Lower values of (3 correspond to larger 
values of ice concentration, as we saw in Fig. 1141 Therefore when total ice concentration increases 
the temperature of the flow tends to increase. This can be explained by the fact that during the 
phase transition in the freezing process there is an emission of heat. Therefore more ice production 
corresponds to more heat in the flow and the average temperature tends to increase. 
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Figure 14: Evolution of ice concentration for different values of (3 for Ra = 10 6 , w r = 100 and 
7 = 10 5 



6.4 Free surface flow 

In subsection l6.ll we were interested in the ice produced in the regions far from the upper layer, 
and therefore we adopted the simple boundary condition C, = for the ice concentration at the 
surface. 

In this subsection we want to study if and how the upper boundary conditions affect the ice 
concentration profile and growth. There is still no general consensus on how to obtain no upward 
ice flux at the air-water interface. Some authors [T4"| and ^Hj) impose Neumann condition 

dCi/dz — h Cl with h c ^ 0, in order to represent the seeding effect due to the mass exchange 
with the atmosphere. It is still not clear, however, how to model the flux h c . Other authors [T7| 
propose to account for the ice crystals escaping across the air-water interface by accumulating 
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Figure 15: Evolution of kinetic energy for different values of j3 for Ra — 10 6 , w r = 100 and 7 
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Figure 16: Mean vertical temperature for different values of /3 for Ra — 10 6 , w r — 100 and 7 
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them in a floating ice slab. Here we consider the situation where no ice sheet has formed at 
the water surface, therefore the water surface is a free surface, and we neglect the seeding effect. 
The suitable boundary conditions for this situation are free slip conditions for velocity and total 
vanishing flux for ice concentration: 



w 



0, 



du 
dz 



dv 

dz 



where S is the effective Schmidt number: -| 



0. 



Sc 



dz 



w r SCi — 



The only difference with subsection 



16.11 is that we impose vanishing total ice concentration flux instead of vanishing convective flux. 

In Fig. 1181 we plot the vertical profile of ice concentration averaged over horizontal directions; 
we also plot the one obtained using the boundary conditions of subsection l6.ll We see that the two 
curves differ significantly near the surface, where the new boundary conditions produce a peak of 
ice concentration. Far from the upper boundary however the two profiles look very similar. In Fig. 
I17l we present a snapshot taken during the linear phase of the large downwelling plumes and of the 
ice concentration. Comparison with Fig. 0] confirms that the main difference is the accumulation 
of ice near the water surface. 




Figure 17: Isosurface of negative temperature fluctuations (-0.2, green) and isosurfaces of large 
ice concentration (.lC max , white). Snapshot taken at a temporal station during the exponential 
growth regime for a fast ice formation regime (Ra = 10 6 , 7 = 10 4 , w r = 150). 



We compare next the growth rate of ice concentration for the new and old boundary conditions: 
in Fig. El we present the growth rate of ice concentration as a function of the rise velocity w r . 
We observe that the two curves are similar. When the rise velocity w r is small, we observe that 
the growth rate using the new boundary conditions is larger. This can be explained by the fact 
that these boundary conditions allow production of ice near the surface. 

We can thus conclude that the new boundary conditions only affect the ice produced at the 
top of the layer, but they do not affect significantly the overall ice production rate and the ice 
concentration profile in the deep layers. 



7 Conclusions 

In this paper we have explored the sensitivity of the ice formation process to a set of non- 
dimensional parameters: the ice melting/freezing source term 7, the ratio between the latter 
and the temperature source term /?, and the ice particle rise velocity w r . Regarding the Rayleigh 
number Ra, we have considered moderately turbulent regimes up to Ra = 10 8 . 
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Figure 18: Vertical profile of ice concentration averaged over the horizontal directions at a temporal 
frame during the exponential growth regime. The two curves differ by the boundary conditions 
adopted at the water-air interface. The parameters are Ra = 10 6 , 7 = 10 4 , w r = 150. 
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Figure 19: Growth rate of ice concentration versus rise velocity for different boundary conditions 
(Ra = 10 6 ,7 = 10 4 ) 
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We have modeled the mixture of ice crystals and water as a two-phase medium composed of 
water and ice particles of fixed diameter. The equations describing the system are the classical 
balance equations for mass, momentum and energy in the ice-water mixture, supplemented by an 
advection-diffusion equation for the ice mass concentration with an additional source term, which 
account for the net ice production during melting or freezing process. In order to account for 
the net heat exchange during the same process, an additional source term has been introduced 
in the temperature equation. The simulation of the complete equations, however, is subject to 
numerical instability due to large round-off errors. The origin of such error is the factor Ra/a 
(a is the coefficient of thermal expansion) in the buoyancy term, which, for values of Ra of 
interest (Ra > 10 s ), is larger than 10 14 . Therefore, we have developed and applied the Boussinesq 
approximation and we have integrated numerically the set of equations making use of a numerical 
code based on second order finite difference. A dynamic LES model has been used in order to 
cope with the subgrid scales that are not considered in the simulations. 

The sensitivity analysis has shown that the ice particle rise velocity and the ice concentration 
source term coefficient 7 significantly affect the frazil ice dynamics. In the large rise velocity 
regime, the ice crystals are rapidly pushed upward and they cannot act as nuclei for ice formation; 
in the opposite regime of small rise velocity they tend to be transported downwards in warmer 
regions where they tend to melt down. The maximum of ice production is obtained in those 
situations where the rise velocity is of the same order of magnitude of the peak rise velocity w. 
We have developed a simple model which captures the trend of the growth rate as a function of the 
relevant parameters. The model has been derived in the approximation where the ice is produced 
far beneath the water surface, however it is found to be valid also for free surface flow. 

The parameter which plays a key role in fixing the concentration value at the statistically 
steady state is the heat exchange source term coefficient (3. The greater is /3, the lower is the 
domain averaged ice concentration. Moreover, the presence of more frazil ice crystals induces a 
rapid fall in the kinetic energy and this energy lowering increases as the saturated ice concentration 
increases. 

The Boussinesq approximation holds only when frazil ice mass concentration is very small, 
i.e. of the same order of the small non-dimensional parameter a = aAT. So, the Boussinesq 
approximation represents well the initial stage of ice formation, but it breaks down when the 
ice concentration becomes significant. Moreover, at significant ice concentrations, the particles 
interaction and the effective viscosity increment due to the presence of the ice particles should be 
taken into account 

Other effects that were not included in the present study and which we postpone for the future 
are: the inclusion of salinity, the occurrence of different classes of ice crystals with different size, 
the effect of different conditions at the top of the layer, for example a horizontal wind blowing 
over the water, or the presence of an ice shelf. 
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